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CONVERGENT-DIVERGENT NOZZLE 

Ching Y. Loh and K.B.M.Q. Zaman 
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Abstract 

At pressure ratios lower than the design value, 
convergent-divergent (C-D) nozzles often un- 
dergo a flow resonance accompanied by the 
emission of acoustic tones. The phenomenon, 
driven by the unsteady shock within the di- 
vergent section of the nozzle, has been stud- 
ied experimentally by Zaman et al[ 1]. In this 
paper, the space-time conservation element so- 
lution element (CE/SE) method [2, 3, 4] is 
employed to numerically investigate the phe- 
nomenon. The computations are performed for 
a given nozzle geometry for several different 
pressure ratios. Sustained ‘limit cycle’ oscil- 
lations are encountered in all cases. The os- 
cillation frequencies, their variation with pres- 
sure ratio including a ‘stage jump’, agree well 
with the experimental results. The unsteady 
flow data confirm that stage 1 of the resonance 
(fundamental) involves a one-quarter standing 
wave while stage 2 (third harmonic) involves 
a three-quarter standing wave within the diver- 
gent section of the nozzle. Details of the shock 
motion, and the flow and near acoustic field, 
are documented for one case each of stages 1 
and 2. 

1 Introduction 

This paper concerns an aeroacoustic resonance often 
encountered with convergent-divergent nozzles when run 
near ‘transonic’ conditions. The resonance is usually 
accompanied by the emission of intense acoustic tones. 
While a casual observer may easily confuse it with the 
well-known ‘screech tone’, it has been shown to be dif- 
ferent in character as well as origin. The frequency 
of the tone increases with increasing plenum pressure. 
The frequency variation may involve a staging behavior, 
i.e., an abrupt jump in frequency. While odd harmonic 
stages take place at lower pressures, the fundamental 


takes place over a wide range of higher pressures. De- 
pending on nozzle geometry, the fundamental has been 
found to persist to pressure ratios as high as 5. The phe- 
nomenon has been identified and studied experimentally 
by Zaman et al [1]. For a discussion of background, tech- 
nological relevance, and pertinent past work from the lit- 
erature, the reader is referred to the cited reference. A 
discussion of numerical works from the literature appar- 
ently capturing the same phenomenon, and the implica- 
tion of those results, will be deferred to the concluding 
remarks in Section 5. 

During the course of the experimental study a numer- 
ical study was initiated to complement the investigation. 
The ‘CE/SE’ method, to be elaborated shortly, was em- 
ployed because of its past success with flows involving 
shocks and acoustic waves [5]. The calculations were 
first performed for a nozzle geometry and pressure ratio 
that corresponded to an experimental test condition. The 
result was encouraging in that the flow not only exhibited 
a quasi-periodicity but the frequency was also in agree- 
ment with the experimental result. Subsequent calcula- 
tions at two other pressures captured the right trend in 
the frequency variation as well as the stage-jump. These 
promising results prompted a further study. This was 
deemed well-justified not only to shed further light on 
the mechanism of the phenomenon but also to gain con- 
fidence in the CFD methodology. 

The ‘Space-Time Conservation-Element and 
Solution-Element (CE/SE) Method’ [2,3,4] was used as 
a numerical platform in the study. As demonstrated in 
previous papers, the CE/SE method is well suited for 
computing waves in compressible shear flows [5] as well 
as vorticity/shock interactions [6], both being pertinent 
in the phenomenon under consideration. Furthermore, 
based on the novel CE/SE non-reflecting boundary 
conditions (NRBC), it is expected that a small near 
field computational domain would be sufficient and 
the simulation could be focused on the region of most 
importance in the resonance. 
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The objective of this paper is to describe the key re- 
sults of the numerical study. The paper is arranged 
as follows. The axisymmetric CE/SE scheme with an 
unstructured-grid is discussed in Section 2 . The initial 
and boundary conditions as well as the CE/SE NRBC 
are discussed in section 3 . The numerical results are pre- 
sented and compared with experimental results in Sec- 
tion 4 , further discussion and summary are given in Sec- 
tion 5 . 

2 The Unstructured Axisymmetric CE/SE 
Euler and Navier-Stokes Solvers 

In this section the CE/SE numerical scheme is summa- 
rized including (a) the Navier-Stokes CE/SE solver, (b) 
the unstructured grid used, and (c) the treatment of the 
source term and the 2 -D axisymmetric approximation of 
LES (large eddy simulation) applied to the jet flow. The 
basic CE/SE principle and details of the Euler schemes 
can be found in the cited original papers [ 2 , 3 , 4 ]. 

2.1 Conservation Form of the Unsteady 
Axisymmetric Navier-Stokes Equations 

In general, the CE/SE method systematically solves a set 
of integral equations derived directly from the physical 
conservation laws, and hence naturally captures shocks 
and other discontinuities in the flow. Both dependent 
variables and their derivatives are solved for simultane- 
ously and, consequently, the flow vorticity can be ob- 
tained without reduction in accuracy. Non-reflecting 
boundary conditions (NRBCs) are also easily imple- 
mented because of the flux-conservation formulation. 

Consider a dimensionless conservation form of the un- 
steady axisymmetric Navier-Stokes equations of a per- 
fect gas. Let p,u,v,p, and 7 be the density, stream wise 
velocity component, radial velocity component, static 
pressure, and constant specific heat ratio, respectively. 
The axisymmetric Navier-Stokes equations then can be 
written in the following vector form: 

Ut + F x + G y = Q, ( 1 ) 

where x,y ^ 0, and t are the stream wise and radial co- 
ordinates and time, respectively. The conservative flow 
variable vector U and the flux vectors in the streamwise 
and radial directions, F and G, are given by: 


(lh\ 


(Fi\ 

(Gi\ 


, F = 

f 2 , G = 

G 2 

wl/ 

7 

\F 4 J 

G 3 

\G 4 J 


with 

Ui = p. U 2 = pu , U 3 = pv f 

Ui =p /( 7 - 1 ) + p{u 2 + v 2 )/ 2 . 


The flux vectors are further split into inviscid and viscous 
fluxes: 

F = F{ — F v , G = Gi — G v , 

where the inviscid fluxes are the same as in the Euler 
equations: 

Fi 1 = U 2 , 

Fi 2 = (7 - 1 ) 1^4 + [(3 - 7 )Ul - (7 - 1 )U%] / 2 lh, 

Fiz = U2U3/UU 

Fi 4 = 7 tW^i - (7 - 1)I>2 \Ul + C/ 3 2 ] /2Ul 
Gn = U 3 , G i2 = UtUz/Ui, 

G a = (7 - 1 )U A + [(3 - 7 )C/ 3 2 - (7 - l)C/f] /2C/,, 
G u = 7C/3C/4/C/1 - (7 - l)U 3 [U 2 + C/ 3 2 ] /2U 2 , 
and the viscous fluxes are: 

F v i = 0 , F v 2 = p{2 u x — §V-V), 

F v 3 = fJ>(v x +%), 

F v 4 — p^hiVjx + (u y + Vx'jv 

7 d U 4 u 2 +v 2 

T^dy^lh 2 

G V 1 — G v 2 — p(v x + 'Uy)-) 

G v 3 = fJb( 2 Vy - |V-V), 

G v 4 = p[2vv y + (u y + v x )u - |(V-V)u+ 

7 d U 4 u 2 +v 2 

J^'dy^Th 2 

where u,v,u x ,u y ,Vx,Vy are respectively the x— and y— 
flow velocity components and their derivatives. They 
can be written in terms of the conservative variables 
I 7 i, 1/27 U 3 and U4. Pr is the Prandtl number, and p the 
viscosity. The velocity divergence 

V-V = u x + v y + v/y. 

The right hand source term Q is the same as in the 
axisymmetric Euler equations [ 7 ]: 

f Ql \ 

Q- ® 2 
Q Qs * 

\q 4 / 

where 


Qi — /y, Q2 — —G 2 lJ 3 /U\y, 

Qz = -U 2 /U iy , Qi = -Gi/y. 

By considering (x,y,t) as coordinates of a three- 
dimensional Euclidean space E s and using Gauss’ diver- 
gence theorem, it follows that Eq. ( 1 ) is equivalent to the 
following integral conservation law: 

H m -dS= [ Q m dV, m = 1 , 2 , 3 , 4 , ( 2 ) 
Js(V) Jv 

where S(V) denotes the surface around a volume V in 
P'3 an d H m = ( F m 7 G m: U m ). 
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2.2 Unstructured Grid for CE/SE 

The CE/SE scheme is naturally adapted to unstructured 
triangular grid. The unstructured version of the CE/SE 
scheme can be briefly described using Fig. 1. Here, 
A ABC is a typical triangular cell centered at O, and 
D,E,F are the centers of the neighboring triangular 
cells where the flow data at the previous time step are 
given. Each triangle center and its three neighboring 
triangle centers form three cylindrical quadrilateral con- 
servation elements or CE s, as shown in Fig. 1. In the 
space-time E 3 space, (2) is applied to the hexagon cylin- 
der ADBECF ( the volume V) that consists of these 3 
quadrilateral CE s. The discrete approximation of (2) is 
then 


/ 

Js(V) 


H, 


■dS = V(Q m )] +1 , (3) 


for m = 1, 2, 3, 4, where H* m = (F^,G* m ,U^). The 
right-hand side of (3), in general, is the volume V times 
the ‘source’ term evaluated at an appropriate gaussian 
quadrature node. Here, the gaussian quadrature node is 
the center of the hexagon ADBECF at the new time 
level. 

In the CE/SE scheme, the above flux conservation 
relation (2) in space-time is the only mechanism that 
transfers information between node points. A conserva- 
tion element CE (here, quadrilateral cylinders) is the fi- 
nite volume to which (2) is applied. Discontinuities are 
allowed to occur in the interior of a conservation ele- 
ment. A solution element SE associated with a grid 
node ( e.g ., D,E,F in Fig. 1) is here a set of inter- 
face planes in E 3 that passes through this node (e.g. 
DAA!D\ DBB’D', EBB'E\ ECC'E', etc . ). 

At time level n, the solution variables U , U x , and U y 
are given at these three nodes. We are to solve U , U x 
and U y at O' at the new time level n + 1 . 

In principle, each of the 3 CE s provides 4 scalar equa- 
tions when (2) is applied to it. There are totally 12 scalar 
equations for the 12 scalar unknowns at O'. The problem 
is solvable. All the unknowns are solved for based on 
these relations. No extrapolations (interpolations) across 
a stencil of cells are needed or allowed. But in reality, 
U at the hexagon center at the new time level n + 1 is 
first evaluated from (3) and then by Taylor expansion, 
U at the center O' of triangle ABC can be obtained. 
Consequently, U x and U y are solved with the necessary 
numerical dissipation added. Details can be found in [4]. 

An important issue is how to accurately calculate the 
surface fluxes of the SE s. For this purpose, within a 
given solution element SE(j,n), where j, n are the node 
index, and time step respectively, the flow variables are 
not only considered continuous but are also approxi- 
mated by linear Taylor expansions: 


U* (x,y,t;j,n) = U] + (U x )](x - Xj ) + 



hexagon cylinder ADBECFA - A’D’B’E’C’F’A’ and 
its 3 CEs, OADB-O’A’D’B’, OBEC-O’B’E’C’, 
OCFA - O’C’F’A’ 


Figure 1 : CE/SE unstructured grid 


(U v )?(y- yj ) + (U t )](t-t n ), (4) 

F*(x,y,t-,j,n ) = F] + (F x )J(x - Xj )+ 

(F y )”(y- yj ) + (F t )](t-t n ), (5) 

G*(x, y , t j, n) = G] + (G x )?(x - Xj )+ 

(G y )?(y- yj ) + (G t )?(t-t n ), (6) 

where j is the node index of D,E or F. The partial 
derivatives of F and G can be related to the correspond- 
ing derivatives of U by using the chain rule, and U t can 
be obtained from (1). Now, the surface flux can be cal- 
culated accurately by first evaluating the flux vectors at 
the geometrical center of the surface through the Taylor 
expansions (4-6). With unstructured grids, the CE/SE 
procedure is simplified and more adapted to complicated 
geometry. Only a single set of mesh points is needed as 
compared to the spatially staggered mesh for structured 
grids and the time-marching is completed in one step 
rather than two. Also, the simple non-reflecting bound- 
ary conditions described previously [5-8] still work 
well with an unstructured grid. More details about the 
unstructured CE/SE method can be found in [4]. The 
weighted a — e CE/SE scheme is used here. 

2.3 Treatment of the Source Term 

The treatment is identical to the one used in [7] and is 
briefly reiterated here. As the source term Q = Q(U) 
itself is a function of the unknown U at the new time 
level, a local iterative procedure is needed to determine 
U. The discretized integral equation (3) reduces to the 
form 

U - Q(U)At = U H7 (7) 
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where Ur is the local homogeneous solution (i.e. the 
solution for Q — 0 locally). Note that Ur only de- 
pends on the solution at the previous time step, i.e. Ur 
is obtained using explicit formulas. A Newton iterative 
procedure to determine U is then 

uii+1) = V (i) _ (|i)-i[^ (C /W) _ Uh]j 
where i is the iteration number and 


$(U) = U - Q(U)At. 


Normally, U at the previous time step is a good initial 
guess U(°) and the procedure takes about 2-3 iterations 
to convergence. 

2.4 Large Eddy Simulation (LES) 

A simplified LES procedure similar to those used in [9] 
is adopted here to account for the strong momentum ex- 
change in the jet shear layer outside the nozzle. In this 
2-D approximation of LES, a simple Smagorinsky’s sub- 
grid scale model is used for the eddy viscosity: 

IH = ( C s A) 2 (2S ij S ij ) 1 / 2 , 


where 


Jij 


i ( 9uj duj 
2 dxj dxi 


A = (AxAy) 1 / 2 , C s — 0.1. Then fi t + ji repalces fi 
in the Navier-Stokes CE/SE solver. 


3 Computational Domain, Initial and 
Boundary Conditions 

As stated earlier, an axisymmetric 2-D CE/SE Navier- 
Stokes solver was used. It should be noted that ‘high 
Reynolds number’, ‘inviscid’ calculations were first per- 
formed for the nozzle internal flow (that is, ji was set to 
zero for the viscous flux terms in the governing equa- 
tions, Section 2). With a uniform flow at the nozzle inlet 
and no-slip condition on the wall, trunction and other er- 
rors in the numerical process imposed a viscous effect 
and there was a boundary layer growth due to this ‘nu- 
merical viscosity’. Physical viscosity was later applied 
to simulate different Reynolds numbers, as discussed at 
the end of Section 4. The 2-D approximation of LES 
(Section 2.4) was applied to the flow outside the nozzle 
in order to properly simulate the free jet. The calculation 
for the nozzle’s internal flow should be deemed more rel- 
evant in the present study. 

The computed domain started at the inlet of the con- 
vergent section of the nozzle as shown in figure 2. The 
domain extended over a rectangular subdomain outside 
the nozzle (21.5 inch in axial direction x 16 inch in ra- 
dial direction; see figure 3). The computation was started 
with the entire flow at rest and with the desired plenum 



Figure 2: Geometry of the C-D nozzle and unstructured 
grid 


pressure applied at the inlet boundary. A no-slip bound- 
ary condition was imposed on the nozzle wall. For the 
outer subdomain, ambient conditions were applied at 
the upstream inflow boundary, non-reflective conditions 
were applied at the upper and downstream boundaries, 
while a mirror-image reflective (symmetry axis) condi- 
tion was applied at the bottom boundary. 

The geometry of the nozzle can be seen in figure 2. 
The inlet diameter is 1.5” . The diameters at the nozzle’s 
throat and exit are respectively D t = 0.3” and D e = 
0.4”. The throat is located 2” downstream of the inlet. 
The convergent part follows a curved contour r = f(x) 
as noted in the figure, while the contour of the divergent 
part is a straight line. The geometry corresponds to Case 
‘3T2' described in [1]. 

In the experiments it has been shown that the origin 
of the resonance is internal to the nozzle (unlike screech- 
tones). Thus, the numerical investigation concentrates on 
the interior flow of the nozzle and the near field. Further- 
more, the flow and the acoustic fields have been found to 
be of the axisymmetric mode in the experiments. Hence, 
the axisymmetric CE/SE code has been deemed suffi- 
cient, at least as a first attempt, in computing the flow. 
The ambient flow around the nozzle is assumed station- 
ary. 

For convenience, a = 1” is chosen as the length scale. 
(However, the throat- to-exit length, L = 0.75”, is used 
to nondimensionalize data in some of the figures to be 
commensurate with the experiment.) The computational 
domain is a circular cylinder of 24.5a long and 16a in 
radius. The unstructured grid is formed by cutting a rect- 
angular cell into 4 triangles ). These rectangular cells are 
non-uniform and their numbers in the x and y directions 
are 200 and 225, respectively. In a typical case, total 
number of triangles is 114,000. The last 10 cells in the 
streamwise direction have exponentially growing x size 
and serve as a buffer zone to ensure no numerical reflec- 
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Flow at the nozzle inlet is at rest: 



16a 


-24.5a 


Figure 3: Geometry of the entire computational domain 
and unstructured grid. 


tion from the outflow boundary. The grid resolution was 
deemed sufficient since a coarser grid produced the same 
results as discussed shortly. The CE/SE scheme is of the 
a-t type [3,4] with a = 0 and t = 0.5. 

3.1 Initial Conditions 

The pressure pi and the speed of sound a\ in the ambient 
flow are used to scale the dependent variables. Initially, 
the flow of the entire domain is set at the ambient condi- 
tions, i.e., 


Pi — 1 — 1-4, Ui = 0, Vi = 0, pi = 1. 


uo = = 0 . 

No artificial forcing is imposed anywhere. As stated ear- 
lier, no slip condition is applied at the nozzle walls. At 
the symmetry axis, i.e. y — 0, reflective boundary con- 
dition is applied. At the top and outflow boundaries, the 
Type I and Type II CE/SE non-reflecting boundary con- 
ditions as described in the next subsection are imposed 
respectively. 

3.3 Non-Reflecting Boundary Conditions 

In the CE/SE scheme, non-reflecting boundary condi- 
tions (NRBC) are constructed so as to allow fluxes from 
the interior domain to a boundary CE smoothly exit to 
the exterior of the domain. There are various variants 
of the non-reflecting boundary condition and in general 
they have proven to be well suited for aeroacoustic prob- 
lems [5-7]. The following are the ones employed in this 
paper. 

For a grid node ( j, n) lying at the outer radius of the 
domain the non-reflective boundary condition (type I) re- 
quires that 

(U X y = (Uy)? = 0, 

while U 'j is kept fixed at the initially given steady bound- 
ary value. At the downstream boundary, where there 
are substantial gradients in the radial direction, the non- 
reflective boundary condition (type II) requires that 


Note that Pi = 1 corresponds to a dimensional pres- 
sure of 14.4 psi in the experiments. Consequently, the 
conservative flow variables and their spatial derivatives 
can be obtained in an easy way. Initially, all spatial 
derivatives are set to zero. 

3.2 Boundary Conditions 

At the inflow boundary outside the nozzle, the conserva- 
tive flow variables and their spatial derivatives are speci- 
fied to be the same as the ambient flow. The inlet of the 
C-D nozzle is connected to the plenum which provides 
a constant pressure p = po to drive the flow. Accord- 
ing to the experimental conditions, 5 different po values 
(4, 8, 10, 15 and 20 psig) are chosen in the computa- 
tion. The actual non-dimensional po values are respec- 
tively 1.277778, 1, 555556, 1.6944444, 2.0416667 and 
2.3888889. 

Following the experimental condition, the temperature 
in the plenum is assumed to be equal to that in the ambi- 
ent. By using the assumptions of constant total enthalpy 
and isentropic flow, it follows that the density pi at the 
nozzle inlet is related to the ambient pressure pi and den- 
sity pi by 


(U x % = 0 , 

while U'j and (J7 y )” are now defined by simple extrap- 
olation from the nearest interior node j 1 , i.e., 

u'j = uyu 2 {Uy) n = {Uy )yu \ 

As will be observed later, these NRBCs are robust 
enough to allow a near field computation without dis- 
turbing or distorting the flow and acoustic fields. 

4 Results 

The calculations were carried out until a ‘limit cycle’ 
oscillation in the flow field was reached. At this state 
the flow property at a given point in the computational 
domain would undergo a quasi-periodic oscillation with 
varying time. The pressure oscillation at a given point 
in the computational domain for the po =1.694 (10 psig) 
case is shown in figure 4. The underlying periodicity 
in the time history should be apparent upon an inspec- 
tion. The power spectrum corresponding to the data of 
figure 4, shown in figure 5, illustrates the periodicity un- 
ambiguously. The spectrum is clearly characterized by a 
peak at a frequency of about 2600 Hz. A large number of 
time steps (410,000 to 740,000) had to be run in order to 
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0.001 0.002 0.003 0.004 0.005 

t (sec) 

Figure 4: Time history of pressure fluctuation at x/L = 
0.667, y/L = 2.66. 

achieve appropriate resolution on the low frequency end 
of the spectrum. 

Also shown in figure 5 is the power spectrum density 
(PSD) calculated at approximately the same location us- 
ing a coarser grid (80,000 cells). The amplitude of the 
peak is somewhat different from that for the fine grid. 
This is not unexpected since the phenomenon is not ex- 
actly periodic and the physical locations are not exactly 
the same. However, it is clear that the frequency of the 
peak is the same. Thus, the fine grid resolution (1 14,000 
cells) may be considered adequate. 

The oscillation could be detected basically everywhere 
in the computational domain (although in regions corre- 
sponding to ‘nodes’ the pressure oscillation might not be 
clear). The p’ -spectrum obtained from data at a different 
point is shown in figure 6, as another example. Except for 
some difference in the higher frequency peaks the over- 
all spectrum and the dominant peak remains unchanged. 
Note that the data in figure 5 represent near-field acous- 
tic pressure (outside the nozzle) whereas those in figure 6 
represent pressure oscillation within the core of the flow. 
The flow oscillation frequencies are determined in a sim- 
ilar manner through spectral analysis for all five oper- 
ating pressures. These frequencies are compared with 
experimental data in figure 7. 

Here, let us first briefly summarize the experimental 
results (Zaman et al. [1]). The frequency data were ob- 
tained by spectral analysis of a microphone signal. The 
data include ‘screech’ tones as well as the ‘transonic 
tones’, as marked in the figure. The latter phenomenon 
is the subject under consideration. Note that there is a 
staging behavior with the transonic tones. Two stages, 
marked in the figure, are detected with the nozzle un- 
der consideration. From an analysis of data from a large 
number of nozzles, it was inferred that stage (1) repre- 



Figure 5: Power spectral density of fluctuating pressure 
at x/L = 0.667, y/L = 2.66. Solid line: fine grid (1 14,000 
cells), dotted line: coarse grid (80000 cells). 



Figure 6: Power spectral density of fluctuating pressure 
within the convergent section at x/L = -0.667, y/L = 0.20. 
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Figure 7: Variation of tone frequency with Mj . Open 
symbols: experimental data (for nozzle 3T2; Zaman et 
al. [1]), Solid symbols: present numerical results. 


sented the fundamental in the resonance while stage (2) 
represented the next odd (3rd) harmonic. It was found 
that the presence of a shock within the upstream reaches 
of the divergent section was a necessary condition for the 
resonance to take place. A physical model for the under- 
lying mechanism was proposed. In short, the unsteady 
shock was thought to act like a vibrating diaphragm and 
resonance took place in a manner similar to that occur- 
ring in the simple (no-flow) acoustic resonance of a con- 
ical section with one end open and the other end closed. 
Thus, the fundamental was expected to involve a one- 
quarter (wavelength) standing wave within the diverg- 
ing section while the next (third) harmonic (stage 2) was 
expected to involve a three-quarter standing wave. Un- 
steady flow measurements did indicate the presence of 
such standing waves. However, the latter observation 
was not on firm ground because of the possibility of 
probe interference effects. 

In figure 7, the numerical results for the resonant fre- 
quencies are shown by the solid symbols. The results 
agree with the experimental data very well. Not only the 
trend of frequency variation within each stage, and the 
stage jump, are captured but also the frequencies are pre- 
dicted quite well. 

Details of the flow field were computed for two pres- 
sures, po =1.278 and 1.694, corresponding to stage (2) 
and stage (1) resonance, respectively. The static pressure 
distributions are shown in figure 8 for the fundamental 
case at po =1.694. The 11 frames span approximately 
2T, T being the period of the oscillation. Thus, approx- 
imately on the 6th frame the period is completed and a 
similar distribution is expected as in the 1st frame. This 



Figure 8: Detailed flow-field pressure data for the funda- 
mental case at 10 psig (stage 1). 


is indeed the case as can be seen. Numerical Schlieren 
pictures (i.e., distributions of |^), corresponding to the 
instants of figure 8, are shown in figure 9. These illus- 
trate the unsteady shock structure and its motion within 
the period. A shock (denoted by the boundary between 
the yellow/red and blue/green regions) can be seen past 
the throat of the nozzle. The shape and structure of the 
shock changes widely over the period. During part of the 
cycle, a clear ‘bow-shaped’ front is seen (frames 1,6). 
During other parts of the cycle a ‘lambda-shock’ is seen 
clearly (frame 2,7). Yet during other parts of the cycle, 
multiple fronts are noted. Corresponding distributions of 
the axial velocity (U) are shown in figure 10. A flow sep- 
aration downstream of the throat of the nozzle can be ob- 
served. The length of the separated flow region changes 
over the period. 

Flow field details for the po = 1.278 (stage 2) case are 
shown in figures 11-13, in a similar manner as in figures 
8-10. Here, an approximate repetition of the flow pattern 
every sixth frame can also be observed. However, it is not 
as clear as in the case of the fundamental. In any case, the 
completion of the period can be inferred, for example, 
by following the yellow patch of high-pressure region to 
the right of the nozzle exit in figure 11. It propagates 
downstream with increasing time step, until in the sixth 
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Figure 9: Detailed numerical Schlieren data for the fun- 
damental case at 10 psig(stage 1). 



Figure 10: Detailed flow-field u- velocity data for the fun- 
damental case at 10 psig (stage 1). 


Figure 1 1 : Detailed flow-field pressure data for the case 
at 4 psig (stage 2). 


frame a pattern similar to that in frame 1 reappears. 

Referring back to the Schlieren pictures for the po 
=1.694 case, a qualitatively similar shock motion may be 
noted in past experimental investigations. The Schlieren 
pictures in figure 14 are from the experiment of C.A. 
Hunter of NASA Langley Research Center for a 2-D 
convergent- divergent nozzle (see Hunter 1998 [10], and 
Zaman et al. [1] for further details). Compare for exam- 
ple, the pictures in figures 14(a), (b) and (c) with frames 
1, 2 and 4 of figure 9. 

The r.m.s. pressure fluctuation amplitudes were com- 
puted from the data on the centerline of the nozzle. The 
results are shown in figures 15(a) and (b) for the cases 
of the 3rd harmonic (stage 2) and fundamental (stage 1), 
respectively. The throat is located at an abscissa value of 
0 while the nozzle exit is at 1. First, in figure 15(b), it 
can be noted that there is a pressure anti-node somewhat 
downstream of the throat while a node exists somewhat 
downstream of the exit. This clearly indicates the pres- 
ence of a one-quarter standing wave within the divergent 
section. Similarly, a three-quarter standing wave is ob- 
served for stage (2) in figure 15(a). The results of fig- 
ure 15 agree with and confirm the presence of the stand- 
ing waves that were conjectured from experimental evi- 
dence, as discussed earlier. 
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Figure 13: Detailed flow-field u- velocity data for the case 
at 4 psig (stage 2). 



Figure 14: Schlieren pictures of the internal shock struc- 
ture for fundamental resonance at Mj = 0.71 for a rectan- 
gular nozzle, from the experiment of Hunter [10]. 


A recent experimental investigation in connection with 
flow metering using Venturi nozzles [11] is worth men- 
tioning here. In this work, flow unsteadiness was ob- 
served at relatively low pressure-ratios. Time-averaged 
data on ‘recovery temperature’ clearly exhibited the pres- 
ence of one-quarter, three-quarter and even five-quarter 
standing waves within the divergent section. (In private 
communication, the first author of [11] confirmed that 
the flows were accompanied by emission of tones. Un- 
fortunately, frequencies were not measured that would 
have allowed a comparison with the correlation equa- 
tions of Ref. 1 and a determination if the unsteadiness 
was indeed the same phenomenon as studied here.) 

In figure 16, the pressure amplitudes are shown for the 
fundamental case at a given instant. The acoustic radia- 
tion pattern at the resonance frequency is captured in this 
plot. It is noteworthy that the pattern is similar to that 
observed with duct acoustic resonance. For example, in 
the work of [12] the pressure fluctuations for a resonating 
cylindrical duct were calculated. While the internal pres- 
sure amplitudes apparently showed standing waves (their 
Fig. 13), the external radiation pattern appeared similar 
to that observed here. 

Finally, the effect of Reynolds number is examined on 
the resonance frequency. In the experiments the reso- 
nance tended to disappear with increasing pressure ra- 
tio (or Mj). Thus, the highest pressure-ratio ( 20 psig 
case) was chosen for this Reynolds number sensitivity 
study. Computations were performed for three additional 
Reynolds numbers, calculated on the basis of acoustic 
speed in the ambient and a length- scale of 1 inch. The 
solid curve, representing ‘high Re’ calculation was ob- 
tained in the same way as described above (by setting 
{i to zero; see Section 3). This involved boundary layer 
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with growth only due to ‘numerical viscosity’. The latter 
depends on grid size and other parameters of the com- 
putational procedure. An estimate of the numerical vis- 
cosity is not straightforward, however, the corresponding 
Reynolds number is thought to be quite high. The lower 
Reynolds number cases were obtained by setting corre- 
spondingly higher values of fi. Note that the nominal 
Reynolds number based on molecular viscosity is about 
600,000. 

The power spectral density of the fluctuating pressure 
is shown in Fig. 17. For brevity, these calculations were 
performed only for limited lengths of time so that the 
resolution in Fig. 17 is not as refined as in Figs. 5 and 
6. The ordinate in Fig. 17 has a log scale and thus the 
higher harmonics appear prominently. It is clear that the 
resonance persists over a wide range of Reynolds num- 
bers and for high and moderate Re’s, their effect on the 
resonance frequency is quite small (about 2% at maxi- 
mum). Only at the lowest Re the oscillation disappears. 
Note that a deceasing Re is equivalent to a thickening 
of the boundary layer. The trend is therefore consistent 
with experimental observation that the resonance is sup- 
pressed upon tripping of the boundary layer. However, 
the exact mechanism of the tripping effect remains un- 
clear [1] and further study, including computational, is 
needed to obtain a full understanding. 

5 Discussion and Summary 

It is important to discuss certain past work that ad- 
dressed the same or similar unsteady phenomenon. In the 
work reported in [13], a simulation of the experiments 
of [14] was performed. The experiment [14] reported a 
periodic unsteady flow through a two-dimensional ‘tran- 
sonic diffuser’ . The floor of the diffuser was flat, the two 
side- walls were parallel and the top wall was convergent- 
divergent. An examination of the frequency data from 
this experiment (see Ref. [1] for full discussion), led 
to the inference that the observed unsteadiness must 
be of same morphology as of the subject phenomenon. 
Of relevance here is the fact that the numerical simu- 
lation of [13] appeared to have captured the flow un- 
steadiness quite well. A two-dimensional, compress- 
ible, Reynolds-averaged Navier-Stokes (RANS) solver 
was used together with a two-equation turbulence model. 
The computation was started with an initial state deter- 
mined from quasi-one-dimensional analysis. After suffi- 
cient time steps the flow settled into an oscillatory pat- 
tern. The period of oscillation for some of the cases 
agreed well with the data of [14]. However, it was found 
that the diffuser flow field and the frquency were very 
sensitive to the location of the downstream boundary. 

Another numerical work [15] brought to the authors’ 
attention also merit a discussion. This work concerned 
flow metering using Venturi nozzles. Numerical simu- 
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Figure 15: Time-averaged (r.m.s.) amplitude of fluctu- 
ating pressure along centerline of nozzle; (a) po = 1.278 
case (stage 2); (b) po = 1.694 case (stage 1). 
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lation was conducted for flows similar to that of [11], 
discussed in connection with Fig. 15. An axisymmetric, 
unsteady, compressible Navier-Stokes solver was used. 
For certain conditions, a quasi-periodic fluctuation was 
observed in the computed flow field. From the data pro- 
vided in [15] it was not possible to determine if the fre- 
quencies followed the correlation equations of [1]. How- 
ever, the symptoms appeared similar and the unsteadi- 
ness seemed likely to be of the same origin as the tran- 
sonic resonance. 

Thus, to the authors’ knowledge, at least two past 
numerical works also likely captured the transonic res- 
onance phenomenon. The apparent ability to capture 
the phenomenon by ‘any’ unsteady code is surprising 
because experimentally it is known to be very sensi- 
tive to operating conditions especially to perturbations 
in the upstream boundary layer (see discussion of Fig. 
17). Except for the similarity that all three simulations 
(present, [13] and [15]) involved axisymmetric or two- 
dimensional codes, different algorithms and procedures 
were followed. The upstream boundary layer in the 
present simulation was ‘laminar’, that in [13] was appar- 
ently turbulent. We move on by noting that many aspects 
of the phenomenon, e.g., the effect of boundary layer 
tripping, have remained far from completely understood. 
Further numerical study has the promise of advancing the 
understanding, and this is planned for the future. Follow- 
ing is a summary of the results presented in this paper. 

The unstructured CE/SE Navier-Stokes/Euler solver is 
applied to the transonic resonance phenomenon. It is 
clear that the essence of the phenomenon is captured by 
the computation. The frequency of the resonance and 
its variation with pressure-ratio, including a stage jump, 
are captured quite well. A significant contribution of the 
present study is the results clearly showing the charac- 
teristic standing waves (Fig. 15). This confirmed that 
the underlying mechanism is similar to that of acoustic 
(no-flow) resonance of a duct having one end open and 
the other closed. This was conjectured in the earlier ex- 
perimental work [1] but the standing waves could not be 
measured with confidence because of probe interference. 

The success in capturing the characteristics of the phe- 
nomenon attests to the validity of the numerical scheme. 
Other advantages of the CE/SE scheme include the ‘ef- 
fortless’ implementation (no special treatment, grid re- 
finement, etc.), the simple but effective NRBC and its 
shock-capturing capability. 

As stated in Section 3, most of the results in this paper 
were obtained by ‘inviscid’ calculations for the internal 
flow. It can be viewed as a Navier-Stokes solution at 
a high Reynolds number with under-resolved boundary 
layer; the computed flow involved a boundary layer due 
to numerical viscosity only. Boundary layer separation 
following the shock can be observed in figures 10 and 13. 



Figure 16: Acoustic radiation for the fundamental case. 
(lOpsig) 


However, one may not expect that the details of the sep- 
arated flow and the separation bubble would be captured 
faithfully by such a procedure. On the other hand, users 
of the CE/SE method have demonstrated its capability to 
capture shock structure relatively faithfully. Thus, one 
may infer that the shock and its unsteadiness, possibly 
due to flow separation, are the primary ingredients of the 
transonic resonance phenomenon and that the exact de- 
tails of the separated boundary layer are here relatively 
unimportant. 


References 

[1] Zaman, K.B.M.Q., Dahl, M.D. and Bencic, T.J. 
“Experimental Investigation of ‘Transonic Reso- 
nance’ with Convergent-Divergent Nozzles”, AIAA 
Paper 2001-0078 (2001). 

[2] Chang, S. C., “The Method of Space-Time Con- 
servation Element and Solution Element — A New 
Approach for Solving the Navier-Stokes and Euler 
Equations,”/. Comput. Phys, vol. 119, pp. 295-324 
(1995). 

[3] Chang, S.-C., Wang, X.-Y. and Chow, C.-Y., “The 
Space-Time Conservation Element and Solution 
Element Method — A New High Resolution and 
Genuinely Multidimensional Paradigm for Solving 
Conservation Laws,” /. Comp. Phys ., vol. 159, pp. 
89-136(1999). 


NASA/TM— 2002-21 1324 


11 




Figure 17: Power spectral density of fluctuating pressure 
ai x/L = .667, y/L = 2.667, showing the influence of 
Reynolds number Re (20 psig case). 

[4] Wang, X.-Y. and Chang S.-C., “ A 2-D Non- 
splitting Unstructured Triangular Mesh Euler 
Solver Based on the Space-Time Conservation Ele- 
ment and Solution Element Method” C.F.D. /., vol. 
8, pp. 309-325 (1999). 

[5] Loh, C. Y., Hultgren, L. S. and Chang S.-C., 
“Computing Waves in Compressible Flow Using 
the Space-Time Conservation Element Solution El- 
ement Method,” AIAA /., vol. 39, pp. 794-801 
(2001); also AIAA Paper 98-0369 (1998). 

[6] Loh, C. Y., Hultgren, L. S., Chang, S.-C. and Jor- 
genson, P. C. E., “Vortex Dynamics Simulation in 
Aeroacoustics by the Space-Time Conservation El- 
ement Solution Element Method,” AIAA Paper 99- 
0359 (1999). 

[7] Loh, C. Y., Hultgren, L. S., and Jorgenson, P. C. E., 
“Near Field Screech Noise Computation for an Un- 
derexpanded Supersonic Jet by the CE/SE Method,” 
AIAA Paper 2001-2252 (2001). 

[8] Chang, S.-C., Himansu, A., Loh, C. Y., Wang, X. Y., 
Yu, S.-T. and Jorgenson, P. C. E. “Robust and Sim- 
ple Non-Reflecting Boundary Conditions for the 
Space-Time Conservation Element and Solution El- 
ement Method”, AIAA Paper 97-2077 (1997). 


[9] Bogey, C., Bailly, C. and Juve, D., “Computation of 
the Sound Radiated by a 3-D Jet Using Large Eddy 
Simulation,” AIAA Paper 2000-2009 (2000). 

[10] Hunter, C.A., “Experimental, Theoretical, and 
Computational Investigation of Separated Nozzle 
Flows”, AIAA Paper 98-3107 (1998). 

[11] Ishibashi, M. and Takamoto, M., “Discharge co- 
efficients of critical nozzles with step near the 
throat and their flow field estimated from recov- 
ery temperature distribution”, FEDSM2001-18035, 
ASME Fluids Engineering Division Summer Meet- 
ing, New Orleans, LA. (2001). 

[12] Hu, F.Q. and Manthey, J.L., “Application of PML 
absorbing boundary conditions to the benchmark 
problems of computational aeroacoustics”, NASA 
CP 3352 , Second CAA Workshop on benchmark 
problems, edited by C. K. W. Tam and J. C. Hardin, 
June 1997. 

[13] Hsieh, T. and Coakley, T.J., ’’Downstream boundary 
effects on frequency of self-excited oscillations in 
transonic diffuser flows,” AIAA Paper 87-0161, 25th 
Aerospace Sciences Meeting, Reno, NV, Jan. 12- 
15, 1987. 

[14] Bogar, T.J., Sajben, M. and Kroutil, J.C., “Charac- 
teristic frequencies of transonic diffuser flow oscil- 
lations,” AIAA /., 21 (9), pp. 1232-1240., 1983. 

[15] von Lavante, E. Zaichcial, A., Nath, B. and Di- 
etrich, H., “Numerical and experimental investiga- 
tion of unsteady effects in critical Venturi nozzles”, 
Flow measurement and instrumentation , vol. 11, 
pp. 257-264, 2000. 


NASA/TM— 2002-21 1324 


12 


REPORT DOCUMENTATION PAGE 

Form Approved 
OMB No. 0704-0188 

Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, 
gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this 
collection of information, including suggestions for reducing this burden, to Washington Headquarters Services, Directorate for Information Operations and Reports, 1215 Jefferson 
Davis Highway, Suite 1204, Arlington, VA 22202-4302, and to the Office of Management and Budget, Paperwork Reduction Project (0704-0188), Washington, DC 20503. 

1. AGENCY USE ONLY ( Leave blank) 

2. REPORT DATE 

January 2002 

3. REPORT TYPE AND DATES COVERED 

Technical Memorandum 

4. TITLE AND SUBTITLE 

5. FUNDING NUMBERS 


Numerical Investigation of “Transonic Resonance” With a 
Convergent-Divergent Nozzle 


6. AUTHOR(S) 


Ching Y. Loh and K.B.M.Q. Zaman 


WU-7 08-90-43-00 


7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 

National Aeronautics and Space Administration 
John H. Glenn Research Center at Lewis Field 
Cleveland, Ohio 44135-3191 


8. PERFORMING ORGANIZATION 
REPORT NUMBER 


E-13134 


9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 

National Aeronautics and Space Administration 
Washington, DC 20546-0001 


10. SPONSORING/MONITORING 
AGENCY REPORT NUMBER 


NASA TM— 2002-21 1324 
AI A A-2002-007 7 


11. SUPPLEMENTARY NOTES 


Prepared for the 40th Aerospace Sciences Meeting and Exhibit sponsored by the American Institute of Aeronautics and 
Astronautics, Reno, Nevada, January 14-17, 2002. Ching Y. Loh, Taitech Inc., 21000 Brookpark Road, Cleveland, Ohio 
44135; and K.B.M.Q. Zaman, NASA Glenn Research Center. Responsible person, K.B.M.Q. Zaman, organization code 
5860, 216-433-5888. 


12a. DISTRIBUTION/AVAILABILITY STATEMENT 

Unclassified - Unlimited 
Subject Categories: 02, 64 and 71 

Available electronically at http :// gltrs. grc .nasa. gov/GLTRS 

This publication is available from the NASA Center for AeroSpace Information, 301-621-0390. 


Distribution: Nonstandard 


12b. DISTRIBUTION CODE 


1 3. ABSTRACT (Maximum 200 words) 

At pressure ratios lower than the design value, convergent-divergent (C-D) nozzles often undergo a flow resonance 
accompanied by the emission of acoustic tones. The phenomenon, driven by the unsteady shock within the divergent 
section of the nozzle, has been studied experimentally by Zaman et al. In this paper, the space-time conservation element 
solution element (CE/SE) method is employed to numerically investigate the phenomenon. The computations are per- 
formed for a given nozzle geometry for several different pressure ratios. Sustained ‘limit cycle’ oscillations are encoun- 
tered in all cases. The oscillation frequencies, their variation with pressure ratio including a ‘stage jump’, agree well with 
the experimental results. The unsteady flow data confirm that stage 1 of the resonance (fundamental) involves a one- 
quarter standing wave while stage 2 (third harmonic) involves a three-quarter standing wave within the divergent section 
of the nozzle. Details of the shock motion, and the flow and near acoustic field, are documented for one case each 
of stages 1 and 2. 


14. SUBJECT TERMS 

Transonic resonance; C-D nozzle; CE/SE method 


15. NUMBER OF PAGES 

18 


16. PRICE CODE 


17. SECURITY CLASSIFICATION 
OF REPORT 

Unclassified 


18. SECURITY CLASSIFICATION 
OF THIS PAGE 

Unclassified 


19. SECURITY CLASSIFICATION 
OF ABSTRACT 

Unclassified 


20. LIMITATION OF ABSTRACT 


NSN 7540-01-280-5500 


Standard Form 298 (Rev. 2-89) 
Prescribed by ANSI Std. Z39-18 
298-102 


